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We describe a numerical scheme for exactly simulating the heat current behavior in a quantum 
harmonic chain with self-consistent reservoirs. Numerically-exact results are compared to classical 
simulations and to the quantum behavior under the linear response approximation. In the classical 
limit or for small temperature biases our results coincide with previous calculations. At large bias 
and for low temperatures the quantum dynamics of the system fundamentally differs from the close- 
to-equilibrium behavior, revealing in particular the effect of thermal rectification for asymmetric 
chains. Since this effect is absent in the classical analog of our model, we conclude that in the 
quantum model studied here thermal rectification is a purely quantum phenomenon, rooted in the 
quantum statistics. 
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I. INTRODUCTION 

Understanding the role of quantum effects in the ther- 
mal conduction properties of interacting systems is a 
challenging task. While in the classical regime molecular 
dynamic simulations provide a flexible tool for including 
anharmonic interactions to all orders 0, [|J , in the quan- 
tum limit treatments are typically limited to particular 
parameter domains @. Among the methods developed 
for tracking the quantum behavior of anharmonic sys- 
tems we recall the non-equilibrium Green's function tech- 
nique which is perturbative in the nonlinear interaction 
strength [J, Q and the master equation approach, which 
is limited to models with weak-system bath couplings and 
to systems with few quantum states [y, |7| . More recently, 
mixed classical-quantum molecular dynamics simulation 
tools were developed, valid at relatively high tempera- 
tures It was also demonstrated that a scheme based 
on the Born-Oppenheimer principle could be constructed 
in the context of thermal conduction, useful for studying 
the dynamics in the off-resonance regime Q. Lastly, ex- 
act quantum simulations of the heat current characteris- 
tics can be performed for simplified models only 10]. 

In this paper we consider the problem of heat transfer 
in the quantum harmonic chain model with each inner 
site connected to a self consistent (SC) reservoir. This 
model has been developed with the motivation to in- 
clude nonlinear behavior in an effective way @, [Tl| - [l3| . 
Specifically, here we would like to gain insight on the role 
of quantum effects at low temperatures on the thermal 
properties of 1-dimensional (ID) chains. For brevity, we 
often refer to this model as the " SC model" . It includes 
a linear chain of N beads connected to N independent 
thermal reservoirs, one at each site, see Fig. [1] While the 
temperatures of the reservoirs attached to the first and 
last particles impose the boundary conditions, the role 
of the inner (self consistent) baths is to provide a simple 
scattering mechanism that might lead to local equilibra- 
tion and to the onset of the (diffusional) Fourier's law of 



heat conduction 0, [l3[ . In practice, the temperature of 
these N — 2 internal baths is determined by demanding 
that in steady-state, on average, there is no net heat flow 
between the chain atoms and these reservoirs. 

The classical version of this model has been proposed 
in Refs. [HI, [III an d recently revisited in [H|,[l4|], demon- 
strating that for long chains Fourier's law is satisfied and 
a linear temperature profile is generated. It has been 
also proved that in the classical regime this model can- 
not support thermal rectification, an asymmetry of the 
current under the exchange of the temperature bias, even 
when some spatial asymmetry is provided [l5l . IToT ] . The 
SC model is also of interest in the context of anharmonic 
lattices [13, EH ■ Overall, it is an example for a hybrid 
model, whose time evolution is dictated both by a Hamil- 
tonian term (deterministic), and by stochastic effects. 




FIG. 1: Scheme of a harmonic chain of TV = 5 beads, where 
the inner particles are connected to SC baths. The wiggly 
lines represent harmonic bonds. The temperatures Ti and 
T5 set the boundary conditions; the temperatures T; (I = 
2, 3, 4) are determined by demanding that the leaking currents 
vanish, Fi — 0. The net heat current across the system is 
given by F\ = — F5. 

The quantum analog of the SC model was studied by 
Visscher and Rich who analyzed the limiting case 
of weak-coupling to the SC reservoirs. More recently, 
the model was revisited by Roy and Dhar [2(J [2l| , who 
demonstrated that under the linear response assumption 
and for asymptotically long chains, Fourier's law holds 
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and a temperature dependent thermal conductivity is re- 
alized. An analytical study of the SC model with alter- 
nate masses has revealed the role of quantum effects at 
low temperatures [22], [23| . More recently, a mathematical 
analysis of the mass-graded SC model in the quantum do- 
main has indicated on the onset of thermal rectification, 
beyond the linear response regime [24j |. 

Focusing on the quantum SC harmonic chain model, 
one should note that exact analytic results are limited to 
the linear response regime [2(| [2l[ . The reason is that be- 
yond this limit, for large temperature biases, the self con- 
sistent condition translates into a set of coupled nonlinear 
equations which seem intractable. Since an analytic so- 
lution is missing, in this paper we suggest a numerical 
scheme for exactly simulating the transport properties of 
this model. The method is useful at low and high tem- 
peratures, in equilibrium, and for far-from equilibrium 
situations. It can be also applied onto three-dimensional 
systems. For simplicity, here we confine ourselves to ID 
models. As an interesting application we perform numer- 
ical simulations on the SC harmonic chain model, incor- 
porating a spatial asymmetry. In accordance with previ- 
ous analytic indications [24], we confirm that the system 
rectifies heat in the quantum regime, beyond the linear 
response limit. 

The paper is organized as follows. In Sec. II we present 
the SC harmonic chain model and describe our numerical 
method. We further explain how to calculate the ther- 
mal properties of the model in the classical limit, and in 
the quantum regime, under the linear-response approxi- 
mation. Section III provides some examples for the heat 
current characteristics in different domains, manifesting 
the onset of the thermal rectifying effect in an asymmet- 
ric setting. Section IV concludes. 



II. MODEL AND METHOD 

We now describe the SC model, introduced in Ref. 
[l2| . The chain includes N atoms, where neighbors are 
connected by harmonic links. Each particle is also bi- 
linearly (position-position) coupled to an independent 
thermal reservoir. The temperatures at the end points 
are set to Ti and T^, establishing the boundary con- 
ditions. In contrast, the temperatures of the internal 
reservoirs T; (I = 2, 3, N — 1) are determined in a self 
consistent manner, by requiring that the net heat cur- 
rent flowing into or from the chain through each contact 
I = 2,3, ...,N — 1 vanishes. For a schematic representa- 
tion see Fig. [T] 

The Hamiltonian of the chain Hs (system), its reser- 
voirs H Bl (baths), and the interaction term V Bl is given 
as the sum of quadratic terms, 

N N 

H = H s + J2H Bl +J2 V Bn (1) 

1=1 1=1 



where 

Hs = -^X^MsXs + -Xg$ s Xs, 

H Bl = -Xg [ M B X Bl + -Xg$ Bl X Bn 

V Bl = X T s V Bl X B . (2) 

Here Ms and M B are real diagonal matrices represent- 
ing the masses of the particles in the chain and the bath 
particles, respectively. The quadratic potential energies 
for the chain and baths are given by the real symmetric 
matrices $s and & B[ , respectively. The term V Bl denotes 
the interaction between the chain and the Ith bath. The 
column vectors Xs and X B[ arc the Hciscnbcrg operators 
of the particle displacements about some equilibrium con- 
figuration. In particular, Xs = {X\, X2, Xjy} with 
Xi as the position operator of the Ith particle of the chain. 
The momentum operators are given by X = M~ 1 P, 
where {Xi,Pi} satisfies the usual commutation relation, 
[Xi,P m ] = ih5i ym . 

Since the system is entirely harmonic, a generalized 
Langevin equation for the system displacements can be 
written [2|, |2fJ, |25[ . This is done by formally solving the 
Heisenberg equations of motion (EOM) for the bath op- 
erators, then plugging them into the EOM of the system 
displacements. In the ohmic limit this results in 

M t X t = -(2Xi - AVi - X l+1 ) - 7l X + m . (3) 

Here Mi is the mass of the Ith particle and the force 
constants are taken as unity. The chain-bath coupling 
strengths are enclosed within the friction coefficients 7;. 
The noise-noise correlations, in frequency domain, satisfy 

The steady-state heat current can be obtained by eval- 
uating (two-point) position-momentum correlation func- 
tions |20|. Specifically, it can be shown that the heat 
current from the Ith reservoir into the chain is given by 

Fl = J2wm cW|[GH] ; , m | 2 - 

m=l J ~°° 7F 

x [/(u;,T;)-/(w,T ro )]. (5) 

Here the matrix G is the inverse of a tridiagonal matrix 
with off-diagonal elements equal to -1 and diagonal ele- 
ments 2- Miuj 2 -z7/w, f(uj,T) = [e UJ ' T -l\- 1 is the Bose- 
Einstein distribution. The temperature profile across the 
system is obtained by demanding that 

Ft=0, l = 2,3,....,N-l. (6) 

This condition translates into a set of N — 2 nonlin- 
ear equations, yielding the inner baths' temperatures T;. 
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Plugging the resulting temperatures inside the expression 
for Fi (or equivalently inside Fn) yields the steady-state 
net heat current flowing across the system, 



J = F,= -F, 



N ■ 



(7) 



An exact-analytic solution of Eq. (J6j) is generally not 
accessible. A mathematical analysis has been carried out 
only perturbatively for a specific model with N = 3 [H| ■ 
In the linear response regime (or in the classical limit) 
these N— 2 equations reduce into a set of linear equations, 
which can be readily solved as we explain below. Here, 
with the motivation to treat quantum systems beyond 
linear response, we develop an iterative numerical scheme 
for acquiring the exact solution of Eq. ([6]), i.e. the set of 
the inner-baths temperatures. 
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FIG. 2: Convergence of the inner reservoirs' temperatures 
with increasing number of iterations. N = 10, Pi = 1, Pn = 5 
and 7n = 0.2 for n — 1, N. The main plot shows the tem- 
perature at site number 5. The inset displays the correspond- 
ing convergence of the net heat current, Fi. 



FIG. 3: (left) The currents Fi, I = 1,2..., TV [see Eq. ©], 
at each site, for a chain with N — 10 beads, obtained after 
applying the iterative procedure p + q — 50 times, /3i = 1, 
Pn = 5, 7; = 0.2. (right) Zooming on the internal currents 
Fi to Fq, zero for perfect SC reservoirs. 



guess for the temperature profile Tj-°\ For example, we 

pick the average temperature = (Ti +T N )/2. These 
values are inserted into the right hand side of Eq. ([8]). 

For each site, we then search for the value Tj which 
yields an equality. We do it by calculating the left hand 
side of Eq. ([8]) over a fine grid of temperatures, searching 
for the temperature T/ 1 ^ which minimizes the difference 

I E™=i Si, m ( T l 1} ) ~ £m=i Si, m (T^)\. This process is 
repeated for each of the inner atoms, to obtain the set of 

corrected temperatures Tj . In the next iteration these 
temperatures are used as the basis supposition, for re- 

(2) 

cciving the subsequent corrected profile T ; . Formally, 
at the kth step we solve the following equation N — 2 
times, for each inner site I, 



A. Quantum case: Exact simulations 

The nonlinear equations ([6]) can be numerically han- 
dled by rearranging the expression for Fi [Eq. ([SJ] as 
follows 

N N 

^ Si, m (Ti) = ^ Si >m (T m ), (8) 

rn—l m— 1 

where 

/oo 
dtvu) 2 \[G(u>)]i im \ 2 —f(ui,T). (9) 

Given T± and TV, our goal is to obtain the temperatures 
of the SC baths I) (/ = 2, ...,N — 1). This can be done by 
following an iterative procedure. We first make an initial 



N N 

J2 ^, m (T/ fe+1) ) = J2 ^, m (T«). (10) 

m— 1 m— 1 

The procedure is repeated until we converge the temper- 
ature profile Ti and the current F\. In other words, we 
maintain their values through iterations. One should fur- 
ther verify that the currents F\, (I = 2, N— 1), flowing 
between the inner sites and the SC reservoirs, are negli- 
gible, as we explain next. 

There are two main sources of error in our procedure: 
(i) The frequency integration in Eq. @ is carried out 
numerically, by discretizing energies between a lower and 
upper cutoffs. Selecting a fine frequency step Au; and a 
large energy cutoff lo c 3> wg, with uj as the chain char- 
acteristic frequency, we have verified that our results are 
robust against Aw and lo c . (ii) Equation (flT))) is solved 
on a discretized temperature grid. It is obvious that for 
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a coarse grid the inner reservoirs' temperatures may sig- 
nificantly deviate from the exact SC values and leakage 
occurs |26|]. To control this error we choose a mesh fine 
enough such that |F;|/|Fi| < 1CT 4 for I = 2,3, ...,N - 1. 
In addition, since for long chains the overall-net energy 
exchange between the SC baths and the system may ac- 
cumulate to large values, we also verify in our simula- 
tions that the incoming and outgoing fluxes, and 
|-Fjv| (equal in principle) differ by less than 0.1%. 

In practice, we found that very delicate grids should 
be adopted for reaching a good accuracy for chains with 
TV > 10. We have therefore developed a two-step pro- 
cedure to improve efficiency. In the first step a rela- 
tively rough grid is constructed, ST = (T x - T N )/200, 
and the iterative procedure [Eq. (fTUl) ] is followed to con- 
vergence, in the sense that the temperature profile stays 
fixed through iterations. However, these temperatures 
still deviate from the optimal (SC) temperatures, and 
significant leakage takes place. We denote by p the num- 
ber of iterations in this part. 

In the second step of our procedure an individual 
mesh is constructed at each site by dividing the sector 
[T ; (p) - ST, T^ + ST] into, say 200 elements. Given these 
individualized grids, we iterate Eq. (fTUl) q more times, 
to converge the temperature profile again. Overall, p + q 
iterations are therefore performed, adopting a two-level 
grid. More generally, one could use a hierarchy of temper- 
ature grids, individually constructed around each particle 
for further improving the accuracy of the results and the 
method efficiency. 

Figures [2] and [3] demonstrate the convergence of our 
scheme, as reflected in three quantities: (i) The tem- 
perature of each internal bath should not vary between 
iterations upon convergence. In conjunction, (ii) the cur- 
rent flowing through the system (F\) should remain fixed, 
(iii) Net exchange of energy between the internal reser- 
voirs and the chain should be rudimentary, relative to 
the current crossing the system. It is important to note 
that one can reach convergence with respect to the first 
two criteria, yet the reservoirs may not act as SC baths 
since the temperature grid selected is too rough. 

For testing our method, we consider a uniform chain 
with N = 10 particles of unit mass. The friction is uni- 
form along the chain, j n — 0.2, n = 1,2, ..,N. Fig. [5] 
presents the temperature at the chain center (T5) as we 
repeatedly solve Eq. (fTUth The small jump after 30 iter- 
ations arises due to the re-definition (and refinement) of 
the temperature grid at this point. The inset shows the 
(concurrent) convergence of the heat current F± — —Fjy. 
We also verify that the reservoirs indeed behave as SC 
baths. Fig. [3] displays the currents in the system, in par- 
ticular the leakage currents Fi^i^n after 50 iterations. 
We confirm that the local leakage is smaller by four or- 
ders of magnitude than the net heat current F\ (right 
panel). This assures us that at the end of the iterative 
procedure the reservoirs serve as SC baths. 



B. Quantum linear-response regime and classical 
calculations 

We outline here the process for obtaining the heat cur- 
rent behavior for the SC model in the quantum linear re- 
sponse regime or in the classical domain. In both cases, 
equation ^ reduces into a set of linear equations. 

In the quantum regime under the linear response ap- 
proximation one assumes that temperature differences 
along the chain are small, 7] — T)_i <C TJ, thus the differ- 
ences of Bose Einstein functions (/; — f m ) in Eq. ([5]) can 
be replaced by the derivative (7] — T m ) x df /dT a with 
T a = (Ti + T/v)/2. This approximation is valid close to 
equilibrium, for \T\ — Tjv| <C Ti,Tjv, or for very long 
chains with small local gradients. Under this approxima- 
tion Eq. ([5]) reduces to 

N r°° lo 

Fi = ]T 7 z7 m / ^ 2 |[G( W )] ; , m | 2 - 

Similarly, in the classical limit the quantum statistics is 
replaced by its high-temperature limit, f(oj,T) ~ T/uj 
and Eq. ([5]) becomes 

N f°° ijj 2 

Fi=J2 Vim / rfw|[GH] /:m | 2 — (T, - T m ). (12) 

m=l J -°° n 

Eq. (TTT]) and Eq. (Til]) are both linear in the SC baths 
temperatures. Therefore, we can organize these equa- 
tions as Fi = J2 m Ci,rn{Ti - T m ), with C Lm containing 
the frequency integration. Demanding that Fi — for 
1=2, ..,N — 1, we get the exact solution [l6[ 

T = A~ x v. (13) 

Here A is a diagonal matrix with N — 2 rows for I = 
2, 3, .., N — 1. Its diagonal elements are Y^m=£i Q,m an d 
the nondiagonal elements are given by — C/ jTO . v is a 
vector defined as vt = Ci.\T\ + C^nTn. The vector T 
includes the sought after inner temperatures T2 to Tjv-i- 
Lastly, given the vector T the current F\ can be readily 
calculated. 



III. RESULTS 
A. Quantum effects in thermal conduction 

We recall that the self consistent reservoirs were intro- 
duced as a tool to include in an effective way anharmonic 
processes. Our objective here is to demonstrate novelty 
in transport mechanisms in the deep quantum domain, 
beyond linear response, as a result of the introduction 
of these SC reservoirs. Furthermore, we compare simu- 
lations using the classical, quantum linear response, and 
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FIG. 4: Heat current as a function of chain size using the 
quantum-exact method (o), quantum linear response formal- 
ism (diagonal), and the classical limit (+). 7; = 0.2 for all 
sites, (a) ft = 0.1, ftv = 0.2, (b) ft = 1, ftv = 1.2, (c) 
ft = 1, ftv = 5, (d) ft = 5, ftv = 20. The dotted lines repre- 
sent the exact quantum behavior when the internal reservoirs 
are detached from the chain. 



quantum-exact treatments, unveiling the importance of 
quantum effects at low temperatures, for systems far from 
equilibrium. In the simulations reported here quantum 
results were obtained using the scheme described in Sec. 
II. A. We refer to these calculations as "quantum exact" 
(QE). We used equation (ITT1) to obtain data in the quan- 
tum domain under the linear response approximation, 
denoted by "quantum linear-response" (QLR). The clas- 
sical (C) behavior was acquired using Eq. (fT2")l . The fol- 
lowing parameters were typically used: unit masses and 
unit force constants [see Eq. ([!])]■ inverse temperatures 
P = 1/T ranging between (3 ~ 0.1 — 20. Within these 
parameters one expects to observe an effective classical 
behavior when j3 < 0.2. We also denote the average tem- 
perature by T a — (Ti + Ijv)/2 and the overall bias by 
AT = Ti ~T N . 

Figure [4] displays the heat current as a function of 
size for a uniform chain, comparing data attained from 
the three different methods: QE, QLR and C. We an- 
alyze the behavior in four cases: (a) at high temper- 
atures corresponding to the classical limit, (b) for in- 
termediate temperatures T a ~ ujs and at small tem- 
perature bias AT < T a , corresponding to the linear- 
response regime (ujs is a characteristic system frequency), 



(c) at low temperatures and for far-from equilibrium sit- 
uations, AT/T Q ~ 1 and, (d) at very low temperatures, 
in the deep quantum regime T a <c ws an d at large bias 
AT/T a ~ 1. 

As expected, in the high temperature regime the three 
methods yield the same value (a). At lower tempera- 
tures, yet adopting a small temperature difference, we 
find that QE and QLR calculations agree, while classical 
simulations overestimate the current (b). At low temper- 
atures and large bias, panel (c) demonstrates that QLR 
calculations overestimate the exact value. In order to 
appreciate the role of the SC reservoirs, the dotted line 
further marks the value of the heat current, which is ob- 
tained within a full quantum calculations while nullifying 
the SC reservoirs. We find that in the absence of these 
reservoirs the current remains fixed. This is indeed the 
expected behavior for harmonic chains with a resonance 
thermal conduction mechanism, applicable at high tem- 
peratures [25j . 

Fig. HJd) exemplifies the heat current behavior at ex- 
tremely low temperatures and for a large temperature 
bias, AT/T a ~ 1. Classical results (not shown) are higher 
by an order of magnitude than quantum data. The fol- 
lowing observations can be made: (i) Within the QLR 
and QE methods, the current demonstrates an enhance- 
ment with size up to N ~ 5, followed by a decay. How- 
ever, in the absence of the SC internal reservoirs (dotted 
line) the current systematically increases with size. This 
behavior could be reasoned as follows: With increasing 
chain length the low frequency modes of a periodic lin- 
ear chain are down-shifted [25|], eventually coming into 
resonance with the populated bath modes. This effect 
leads to the enhancement of the heat current with size, 
the behavior indeed detected in the absence of the SC 
baths [25|. However, when the internal SC reservoirs are 
added, scattering mechanisms are responsible for an ef- 
fective diffusional motion, resulting in the decay of the 
current with size Q • The combination of these two trends 
produces the turnover of the current around N ~ 5. (ii) 
In the deep quantum regime the QE current is higher that 
the QLR result, in contrast to the behavior observed in 
Fig. H{c). As we show in Fig. [5jd) below, QLR calcu- 
lations underestimate the factual temperature profile, in 
this case by about 30%. This fact can explain the sim- 
ilar deviation in the current. On the other hand, QLR 
calculations may produce a more significant temperature 
gradient within the chain, see e.g. Fig. EJc), resulting in 
currents larger than the QE data. As a result of these two 
counteracting factors, it is not trivial to predict ad-hoc 
whether the accurate quantum result is above or below 
the QLR limit. 

We now display the temperature profile for a chain 
with iV=10 particles. We calculate it using the three 
methods, QE, QLR and C, in the four parameter do- 
mains mentioned above. Fig. [5] shows that the three cal- 
culations generally agree at high temperatures and for 
AT/T Q <C 1. In contrast, at very low temperatures and 
for AT/T a ~ 1 the QE results deviate from the QLR and 
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FIG. 5: Temperature profile of the SC reservoirs at site n, 
N = 10, using the exact quantum method (o), quantum linear 
response formalism (diagonal), and the classical expression 
(+), ji = 0.2 for all sites, (a) ft = 0.1, ftv = 0.2, (b) ft = 1, 
ftv = 1.2, (c) ft = 1, ftv = 5, (d) ft = 5, ftv = 20. 



the classical data in a profound way. Specifically, for this 
(symmetric) setup QLR and classical calculations provide 
reservoirs temperatures which symmetrically vary around 
the average temperature T a = (7\ + T/v)/2 [lf|. In con- 
trast, QE calculations reveal that the chain temperature 
is actually higher than this temperature. This shifted 
profile stems from the nonlinear Bose-Einstein distribu- 
tion function characterizing the reservoirs statistics. It 
is also of interest to note that classical and QLR simula- 
tions generally predict an internal temperature gradient 
which is higher than the QE value. Specifically, with 
classical simulations we compute a local gradient which 
is larger by an order of magnitude from the QLR and the 
QE behavior [Fig. Eel)]. This considerable disagreement 
reflects itself in the classical heat current, in this case 
high by one order of magnitude compared to the exact 
results. 



B. Application: Quantum thermal rectification 

Thermal rectification, an asymmetry of the heat cur- 
rent for forward and reversed temperature gradients, has 
been extensively analyzed in the last decade [13, . In a 
desirable rectifier the system behaves as an excellent heat 
conductor in one direction of the temperature bias, while 




FIG. 6: A scheme of the mass-graded harmonic chain with 
SC baths for N = 4 beads. 



for the opposite direction it effectively acts as an insula- 
tor. It is agreed that junctions incorporating anharmonic 
interactions with some sort of spatial asymmetry should 
exhibit this effect [13] ■ Since the SC model includes, in an 
effective way, anharmonic interactions through the action 
of the SC reservoirs, it is of interest to explore whether 
this model could demonstrate the rectifying effect when 
some spatial asymmetry is incorporated. This question is 
of particular interest since neither the classical SC model 
nor the QLR SC case can show thermal rectification, even 
when asymmetry is introduced (lR liH f23j ] . In a recent 
paper analytical arguments were put forward, indicat- 
ing that quantum SC systems should rectify heat [24j . 
In what follows we demonstrate that this effect indeed 
exists in asymmetric quantum harmonic systems with 
SC baths. We incorporate asymmetry either by using 
a mass-graded chain, see Fig. [BJ or by connecting the 
chain asymmetrically to the reservoirs at the boundaries. 

Fig. [7] displays the absolute value of the heat current 
for forward (J+) and reversed (J_) temperature biases, 
studying a mass-graded system with Mi = 0.2 and Mi — 
Mi + 0.2 x (I - 1), using the QE method. While at high 
temperatures and for AT <C T a the effect is negligible 
and </+ ~ J_ , at low temperatures and for large bias J+ 
and J_ evidently deviate, with the current being larger in 
the direction of increasing masses. We also find that the 
rectification ratio |J+/J_| is increasing with chain size, 
an observation which can be reasoned by the growing 
mass difference along the chain. 

We note that in different experimental and theoreti- 
cal studies, e.g., Refs. |28l431| . the opposite tendency 
has been reported, and the preferred direction of heat 
transfer occurs from heavy to light atoms. It is clear 
that the preferred direction depends on the details of the 
system studied [29| - |32T |. For example, Ref. [3(| reports 
on molecular dynamics simulations of heat conduction 
in mass-graded chains assuming anharmonic interactions 
between particles. The preferred heat transport direc- 
tion in that model (heavy to light) is attributed to the 
vibrational coupling between low and high modes, aris- 
ing due to anharmonicity in the system. The dynamics 
is further interpreted in terms of the overlap between the 
power spectra at the chain ends. In contrast, simula- 
tions of thermal conduction in mass-graded nanotubes 
showed the opposite trend, explained by the transfer of 
vibrational e nerg y from the transverse to the longitudinal 
direction [33l[34|. 
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In order to better understand the mechanism of ther- 
mal rectification in our model we display the temperature 
profile for the rectifying system in Fig. [8] At high tem- 
peratures (a) there is a reflection symmetry with respect 
to the average temperature, when the temperature bias 
is reversed [16I ]. In contrast, in the quantum domain be- 
yond linear response an asymmetry is discovered (b) : The 
temperature gradient at the chain center is larger when 
the light masses are in contact with the hot bath, than 
the gradient generated in the reversed case. In particular, 
VT 0.0230 (-0.0185) when the heat flows in the direc- 
tion of increasing (decreasing) masses weight. The ratio 
between these gradients indeed fits |J + /J_| for N = 7, 
see Fig. [Tfc). We emphasize that arguments based on 
the power spectra overlap [10, HH cannot be put on for 
the harmonic SC model since it does not include a phys- 
ical mechanism for coupling different vibrational modes. 
This can be seen in Eq. ([3]): The heat current is retrieved 
by integrating over separate contributions, summing dif- 
ferent frequency components. 

One could also generate a spatial asymmetry in a mass- 
uniform SC model by coupling it unequally to the two 
ends, using 71 ^ y^. Since this is a contact asymmetry, 
one would generally expect its effect to diminish with 
size. Fig. [9] still shows a small enhancement of the recti- 
fication strength with TV, probably due to the increased 
importance of the scattering mechanisms (mimicking an- 
harmonicity) with size. 

Concluding this section, an interesting outcome of our 
study is the confirmation that the quantum harmonic 
chain with SC baths acts as a pure quantum thermal 
rectifier, since its classical analog, or the quantum model 
in the linear response regime, can not demonstrate this 
effect [13] . This behavior is attributed to the introduction 
of the SC reservoirs, whose statistics at low temperatures 
is a nonlinear function of the local temperature, unlike 
the high-temperature or linear response behavior. 



IV. SUMMARY 

We developed a numerical method for acquiring the 
heat current and the temperature profile in the quan- 
tum harmonic chain model with SC reservoirs, beyond 
the linear response approximation. While the technique 
is generally valid for 3D models, we applied it here on ID 
linear chains. At low temperatures and for large temper- 
ature biases we found that the exact quantum results sig- 
nificantly deviate from the QLR and classical behavior. 
As an application, we explored the thermal rectification 
effect in asymmetric systems, either by introducing mass 
asymmetry or by imposing a contact asymmetry. In both 
cases we concluded that quantum statistics is responsible 
for the onset of the nonlinear rectifying effect. 

Our method could be generalized in two nontrivial 
ways. First, the scheme could be feasibly extended to 
describe non-ohmic reservoirs, for studying the role of 
memory effects on thermal transport. This could be done 
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FIG. 7: Thermal rectification in a mass graded system: The 
magnitude of the heat current as a function of chain size for 
forward temperature bias Ti > Tjv (o) and for the backward 
direction, Tm > Ti (square), 7„ = 0.2 for all sites, Mi = 0.2, 
M n = Mi + (n — 1) x 0.2. In the forward direction we used 
(a) ft = 0.1, ftv = 0.2, (b) ft = 1, ftv = 1-2, and (c) ft = 1, 
ftv = 5. The opposite polarities were used in each case to 
generate J_. 




FIG. 8: Temperature profile for a mass-graded N — 7 chain. 
In the forward direction (o) we used (a) ft = 0.1, ftv = 0.2 
and (b) ft = 1, ftr = 5. The opposite polarities were used to 
generate the reversed profile (square). Other parameters are 
the same as in Fig. (7) The inset zooms on the chain center. 
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FIG. 9: Thermal rectification in systems with contact asym- 
metry: The magnitude of the heat current as a function of 
chain size for forward temperate bias T\ > Tjv (o) and for the 
backward direction, Tjv > Ti (square) 71 = 0.1, 7jv = 0.4, 
y n -£ 1:N = 0.2. /3i = 1, Pn = 5, and the the reversed contact 
symmetry, (a) Rectification ratio \J+/J-\. (b) The forward 
and backward currents at high temperature /3i = 0.1 and 
13 n = 0.2, and the reversed setup, with no significant rectifi- 
cation observed. 

by introducing frequency dependent friction terms in Eq. 
([5]). Since the method is fully numerical, one can easily 
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